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Abstract - Molecular dynamics simulation of a generic polymer model is 
applied to study melts of polymers with different types of intrinsic stiffness. 
Important static observables of the single chain such as gyration radius or per- 
sistence length are determined. Additionally we investigate the overall static 
melt structure including pair correlation function, structure function and orien- 
tational correlation function. 



1 Introduction 

The detailed structure of bulk amorphous polymers is a topic of scientific inter- 
est because it is necessary for the microscopic understanding of their properties. 
However, because of the amorphous nature of polymer melts and glasses, struc- 
tural information is difficult to obtain experimentally. 

In particular, it is interesting to know how amorphous a polymer melt is on 
a local scale, i.e. how much residual order is left on a local scale, how far local 
order extends before it disappears into the long-range disorder of amorphous 
systems |Q, and how the local order depends on the molecular architecture. 
The interest has recently been revived by solid-state NMR studies of Graf et al. 
[^| , from which it was inferred that a melt of polybutadiene is far more ordered 
than hitherto expected. The alignment of polymer chains is restricted to a local 
scale, there is no sign of nematic ordering. 

In order to get a better understanding of local packing and ordering effects 
computer simulations are very helpful, because the system is precisely known 
and because one has access to all data including positions and velocities of all 
particles at all times. Atomistic simulations may be useful in order to get an 
understanding of a specific system whereas simplified models yield the properties 
of generic polymer melts. Additionally, they need much less simulation time 
which allows to tackle relatively big systems for long times || f|. Therefore, 
a simple bead-spring model may be a good starting point to elaborate generic 
packing effects. 

There was some work done for semiflexible chains, both by Monte Carlo 
and molecular dynamics, mostly to study liquid crystals || ^, |7|, |) or focusing 
on confined systems |)], |T^|. The influence of chain stiffness on the dynamic 
structure factors of polymer melts was investigated also by analytical theory by 
Harnau et. al JIf[ who found major discrepancies to the fully flexible system 
for large scattering vectors, i.e. on short distances. 

In a previous article JllJ, we showed that there is considerable local chain 
alignment even in melts of fully flexible chains (persistence length: I monomer 



diameter) . This persistence length originates from excluded volume interaction. 
If there was no interaction at all (except for connectivity) the persistence length 
would be zero (e.g. polycatenans). In the present contribution, this model 
is extended to include some more information about the chemical architecture 
of the polymer. We firstly introduce bending potentials of different strength, 
in order to study the effect of semiflexibility of single chain structure as well 
as on the mutual local orientation of neighboring chains. Secondly, we study 
models with alternating stiffness in an attempt to mimic simplistically polymers 
with rigid subunits connected by more flexible links, like polybutadiene and 
polyisoprene with their alternating single and double bonds which are currently 
under investigation experimentally ^ . 



2 Simulated System 

We performed molecular dynamics simulations (for details of the parallel pro- 
gram POLY, see ref. Q) of melts of polymer chains at a density p* = 0.85 
and temperature T* = 1 at a timestep 5t* — 0.01 using a truncated and 
shifted Lennard- Jones potential (Weeks-Chandler- Anderson potential) for the 
excluded- volume interaction between all beads. Lennard Jones reduced units 
are used throughout this paper where the mass m, the potential well depth e 
and the radius of the potential minimum a define the unit system. 

+ e, r < r cuto ff = V2a (1) 

and a finitely extendable non-linear elastic (FENE) potential 

V FEN E{r) = ?~ In (\ - ^) , r<R = 1.5a, a = 30 (2) 

for the connection of neighboring beads. Additionally, a bond angle potential 

V angle =x(l- T -^^- \ (3) 
V ri,i+i / 

is used. This model system was already widely studied both for flexible [|[ [H| 
and for semiflexible or liquid crystalline polymer systems [l0| . 

To first approximation, there is -f- — where l p is the persistence length 
(see section 3)and lb the bond length. In our units, the numerical values for x and 
l p therefore coincide. This potential is applied to every bead, to every 2nd bead, 
every 3rd bead etc. The latter is a useful model for polymers with alternating 
stiffness such as single-bond, double-bond sequences or for copolymers with 
different persistence lengths of the constituents. In the following, we refer to 
a system with angular potential strength x and a (topological) distance of y 
monomers between two applications of the bond potential as x-y system. In 
this sense, a fully flexible chain is referred as 0-1 chain. A 5-2 chain for example 
has x = 5 applied to every second bond angle. 



V LJ (r) =4e 
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Figure 1: a) Reorientation of the end-end vector for 5-1 chains (50 monomers), 
b) Reorientation of segments of length d (d + 1 monomers) for 5-1 chains of 
length 200. 



All simulated systems contained 500 chains of 50, 100 or 200 monomers each, 
so the overall number of particles was between 25.000 and 100.000 in a cubic 
periodic box. 

The short-chain systems (50 monomers) could be observed until the auto- 
correlation function of the end-to-end vector R s was decayed. Figure |l|a shows 
the reorientation in the case of the 5-1 system which has the longest relax- 
ation time. The second Legendre polynomial P2(z) — (3z 2 — l)/2 is used for 
consistency with analyses further below. 

At this time, approximately, the mean square displacement of a single monomer 
begins to coincide with that of the center of mass. For the longer chains, we 
first waited until the radius of gyration and the end-to-end distance did not 
change systematically any more but fluctuated only around their mean values. 
The loss of local orientation of shorter chain segments is shown in figure [l]b. 
One sees that there are two regimes. On short time scales (t* < 5000), there 
is a fast decay due to local processes. On long time scales, however, there is a 
long tail which is determined by the overall motion of the whole chain. For the 
investigations in the following this overall motion is not important. 

All systems were simulated at least for t* = 20000. We trust that the static 
chain properties were well equilibrated, because the overall properties like gy- 
ration radius settled and at least local orientation decayed. Moreover, the error 
estimation for R s and R s was performed according to a binning analysis 
The correlation times for the observed properties resulting from this analysis 
were also exceeded substantially. In the case of 5-1 with 200 monomers, which 
has the longest equilibration times, this "binning time" is about tt = 8000. 
Such systems were then simulated for t* = 40000 to 80000. 
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3 Chain Structure in the Melt 



In this section, we investigate the effect of the melt environment on single chains. 
The presence of the other chains screens out the excluded volume interaction 
and the chain statistics of a self-avoiding walk appropriate for chains in good 
solvent turns into a simple random walk This is, for example, evident in the 
single chain structure functions (see figure |^ in section In the semiflexible 
case, one expects that, at large scales, the Gaussian statistics (random walk) 
is fulfilled, whereas on short scales the local stiffness is relevant. Two concepts 
can be used for analysis: One is the idea of a Kuhn length Ik which is defined 
via 

This assumes that the melt consists of "blobs" of length Ik which contain inside 
all the local information which is not relevant on the long scales. 

The second idea is the persistence length l p which derives from the worm- 
like chain model fiL [l5| . It corresponds to the decay length of the correlation of 
bond orientations (the tangent vector) along the chain. 

(cos a(s)) — (u(s)u(O)), s : monomer index (5) 

which can be shown to decay exponentially in this model 

(cosa(s)) = e - s/i "V (6) 

Since we do not only apply the bond angle potential to every bond, but also 
investigated systems with alternating stiff and flexible bonds (e.g. 5-2), the 
persistence lengths of these systems are not a priori known (at least the x-y 
case y ^ 1). In order to determine l p the bond correlation function (eq. ^ was 
determined, in 100 configurations after the equilibration and the initial decay 
was fitted with an exponential e~ l / lp (see table Q). If the bending potential was 
applied to every monomer the decay was well approximated by an exponential 
and the decay length l p was not too far from the expected value x from the bond 
angle potential. This is in agreement with Monte Carlo results for stronger stiff- 
ness ||. In the case of alternating stiffness, minor deviations from exponential 
decay were observed (see figure ^|a). The error in the bond correlation is about 
0.03. Hence, the systems with very short persistence lengths, were difficult to 
determine because only very few points for fitting the decay were available and, 
therefore, the resulting error bars are not negligible. However, in all cases a fit 
over more than one order of magnitude was possible. 

It is not clear if for the x-2 case the bond correlation has to follow an ex- 
ponential law. However, we found this always to be the case. From a simple 
argument, the effective persistence length in the case of persistence lengths l p % 
and l p 2 for alternating angles is 



11/1 1 

lp 2 \ Ipi lp2 



(7) 
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Table 1: Radius of gyration, end-to-end distance (in LJ units), persistence 
length (in monomers) and Kuhn segment length (in L J units) . A x-y system has 
a stiffening potential of strength xksT applied every y monomers. The errors 
are determined via a binning analysis for error estimation Jl4| . 



because 

e -2l/l p _ e -l/lpi e -l/l P 2 _ ^ 

This is exactly true for all points with even monomer distances in the bond 
correlation function. This result may be generalized to a repetitive sequence of 
n different bond angle potentials 



l p n ^ l PJ 

A more elaborate calculation in the framework of a generalized wormlike chain 
model with varying stiffness yields the same result. 

The persistence length of the fully flexible model is found to be exactly 
one monomer distance which is on average h = 0.97. The bond correlation 
functions (of inner monomers) show in the very beginning a decay with a per- 
sistence length which is close to the expected value (on the length scale of about 
5 monomers). This suggests that the very local orientation correlation is de- 
termined by the "true" potential strength whereas on longer scales finite size 
or many chain effects contribute considerably. This is especially reflected in 
the persistence length values for the 5-1 chains (see figure ||b), where the bond 
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monomers s monomers s 

Figure 2: Bond correlation functions: a) Comparison of systems with uniform 
and alternating stiffness. The lines indicate linear least square fits, whose slopes 
define an effective persistence length, b) Different lengths in the 5-f system. 



correlation function for the shorter chain shows at distances s > 5 substantial 
differences to the longer chain. They may be attributed to finite chain length 
effects. All bond correlation functions were determined starting from the inner- 
most monomers in order to avoid end effects as much as possible. 

Also the end-to-end distances and radii of gyration of the corresponding 
chains were calculated. They are also presented in table [I]. Upon increasing x, 
the x-1 systems stretch the chains considerably (see figure ||). A much larger 
bending force constant is needed for the x-2 systems than for x-1 if one wants 
the same R g . Therefore the systems x-2 behave more like fully flexible chains 
with a bigger monomer which is most strongly seen in the persistence length. 
Note that even in the f 00-2 case where an almost rigid and a fully flexible bond 
alternate, the chain stretching is not as strong as in the 2-f case. So there is a 
fundamental difference between these two scenarios. 

At least in the x-1 cases we find 2l p « Ik as expected from the wormlike chain 
model |Q . The relation R 2 w 6R 2 for the Gaussian chain is well fulfilled in most 
of our cases. The larger deviations, e.g. in the 5-1 case with 50 monomers, may 
be attributed to finite chain length effects. Therefore, we do not see substantial 
deviations from Gaussian behavior. 



4 Melt structure 

Local orientation of neighboring chains may be measured by the spatial orien- 
tation correlation function OCF. To this end, we define unit vectors between 
adjacent monomers u = -jr— jr— j- . The scalar products between two such unit 
vectors describe the angle between chain tangent vectors 

COS a{r) = VLchainl ■ VLchain2- (10) 



6 




012345012345 



Figure 3: a) Radius of gyration and b) End-to-end distance for chains of 
length 50 depending on bending strength for x-1 and x-2 systems (uniform and 
alternating stiffness, rcsp.). 



The distance r denotes the distance between the centers of mass of the respective 
chain segments. In order to compare better to NMR experiments, as well as 
to avoid the distinction between head and tail of the chain, we use the second 
Legendre polynomial P%(r) — i(3cos 2 a(r) — 1). 

Figure |]a shows inter-chain orientation correlation functions of different sys- 
tems. The first minimum (r < 1) is close to P2 = — \ which would indicate a 
perfect perpendicular ordering. Two chains which come so close can only pack 
perpendicular because of the excluded volume interaction. The radial distribu- 
tion function (RDF, see below) shows that there are very few such contacts. 
The first peak (r w 1.2) shows a preferred parallel alignment at the distance 
of the first neighbor. A second parallel peak follows at r 2. The intervening 
minima (r ~ 1.6) get weaker for stronger orientation which indicates a stronger 
local parallel ordering. The OCF decays to zero with r because the system is 
globally isotropic, not nematic. 

The local ordering is only slightly different for the systems 0-1, 2-1, and 5-2, 
whereas the 5-1 chain shows a more pronounced local parallel orientation. For 
the 5-1 chains, there is residual parallel ordering even at the intermediate min- 
imum (r w 1.6), where the other systems show some perpendicular ordering. 
Except for the very few direct contacts, there is parallel orientation between 
neighboring chains. This ordering is visible up to about three monomer diam- 
eters. The more flexible systems (0-1, 2-1, 5-2) show qualitatively a similar 
ordering, but it is less pronounced and there is a intermediate preferred per- 
pendicular orientation at the distance of about the first minimum in the radial 
distribution function (r f» 1.6). 

The orientation depends only weakly on the chain length (figure ||d). This 
demonstrates that the effect is strictly local. This holds even though the global 
dynamics of the chains of different lengths is strongly different: The chains of 
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Figure 4: a) Inter-chain orientation correlation functions for chains with 50 
monomers, b) Inter-chain orientation correlation function for different chain 
lengths for the 5-2 system. In order to better distinguish between the lines a 
running average (over 15 points, Sx ~ 0.05) was performed. 



length 50 are not yet entangled whereas the longer chains are already influenced 
strongly by entanglements (entanglement length in the 0-1 case ss 60 monomers 

©■ ' 

Also orientation correlation functions of longer chain segments are investi- 
gated. In this case, not only vectors connecting nearest neighbors but vectors 
connecting next-to-nearest neighbor beads or beads farther apart are taken into 
account (see figure ||) . 

It is clear (figure ^a) that the effect of local parallel chain orientation is 
not restricted to segments of 2 monomers only. It persists when larger chain 
fragments are analyzed. On the other hand, the degree of ordering decreases 
with the segment size considered. Figure [|d shows again the more pronounced 
local ordering in the 5-1 case compared to the more flexible chains. The 2-1 
and the 13-2 systems coincide. Their persistence lengths are quite similar, and 
for the bigger segment sizes the exact local realization of this persistence length 
seems to average out. 

Our results are also in qualitative agreement with an early lattice Monte 
Carlo investigation of shorter chains |^]]. Lattice models, however, are biased in 
favor of orientation correlation. 

The inter-chain radial distribution function g(r) (RDF) on large scales does 
not change much in the case of added stiffness (figure 0a). However, there are 
some differences on very local scales. Both the second and third neighbor peaks 
are farther apart for stronger stiffness. Furthermore, the minimum between 
the first and second neighbor shell is not as pronounced as in the more flexible 
cases. The local stretching allows a closer approach of chains. This leads to a 
reduction of the expected correlation hole. In fully flexible systems the number 
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Figure 5: Definition of unit vectors of bigger segments along the chain. Centers 
of masses are the midpoints of the arrows. 




Figure 6: Spatial orientation correlation functions of segments of length d: a) 
Different segment lengths in a 5-1 system with 50 monomers, b) Different 
systems, d = 2. We applied a running average in order to be able to see 
differences of the curves. 
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Figure 7: Inter-chain monomer- monomer radial distribution functions: a) Dif- 
ferent stiffnesses for length 50. b) Different lengths for 5-1 system. In order to 
better distinguish between the curves, a running average was applied. 



of neighbors of one monomer being on the same chain increase with increasing 
chain length. No effect of chain length is seen here (figure Qd) which again 
reflects the strict locality of the structure formation in the melt. 

5 Structure functions 

The structure of single chains and of the overall melt may be additionally char- 
acterized by static structure functions. Figure || shows the single-chain and melt 
structure functions of our systems. The (isotropically averaged) melt structure 
function is defined as 

, N c n b 

Smelt(k) = -(I ^exp^fcr™)! 2 ) = S SC (k)S mter (k) (12) 
m=l j=l 

where Ssc denotes the single chain structure function 

1 N c n b 

ssdk) = ^E<iE ex p(^ m )i 2 >- ( 13 ) 

m—l j—1 

The first sums run over all chains (Nc- number of chains, m: chain index), 
the second along the chains (n^ = N /Nc'- number of beads along the chain, j: 
monomer index along the chain). 

In the limit for k — *■ the chain structure is no more visible but we just see a 
massive object which is related to the first plateau in Ssc- The next "scaling" 
regime is connected to the fractal nature of the chains. The self similarity 
yields a decay with kr~ , where d — — is the fractal dimension of the chain. 
For a Gaussian chain we have v = i. Stretching leads to a smaller fractal 
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Figure 8: Structure functions for chains of 50 monomers: a) Single-chain 
structure function, b) Melt structure function. 



dimension resulting in a less steep decay. In the large-fe range we deal with 
structure of the size of one monomer. A bead spring model has no structure 
on a shorter scale. The melt structure function is the Fourier transform of the 
density-density correlation function. It shows therefore some additional peaks 
which correspond to peaks in the RDF. Hence, it contains not only information 
about the single-chain structure but also about the overall structure of the whole 
system. 

The single-chain structure functions Ssc of the x-1 systems look very sim- 
ilar, whereas there are strong differences between the 3-1 and the 3-2 system, 
the latter behaving very much like a fully flexible system. The crossover to 
the scaling regime is shifted to smaller fc-vectors for stiffer chains. There is 
additionally a crossover (k — 0.8) between two regimes in the decay which 
means that there is different fractal chain structure on different length scales. 
At larger scales (small-fc regime) the slope does not differ much from the fully 
flexible case, whereas on intermediate scales larger deviations occur, which in- 
dicate local chain stretching. The 3-2 system, however, is close to the fully 
flexible (Gaussian) system. But its slope is not as steep, which hints at a slight 
stretching of the chains compared to the Gaussian chain. This minor difference 
between the alternating stiffness and the fully flexible chains supports our ear- 
lier suspicion that the alternating chains behave like renormalized bead-spring 
chains with larger effective monomers. On the other hand, the stiffened chains 
with true semiflexibility are strongly different on intermediate scales. 

The melt structure functions S me u{k) differ also between the x-1 case and 
the x-2 case. The latter is again very similar to the 0-1 system. The slopes 
in the regimes around k ~ 1 are clearly different, whereas the fine structure 
revealing peaks connected to the neighboring shells is quite similar. The exact 
positions of these peaks differ, however, which shows again that the distance to 
the nearest neighbors is slightly altered with stiffness. Figure ||b shows that all 
differences are on local scales (k > 0.1). The static structure functions coincide 
for k — ► 0, where the overall structure on large scales gets important. So the 
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melt structure function shows that all systems behave similar on large scales but 
the systems with homogeneous and alternating stiffness differ on local scales. 

6 Conclusions 

The static structure of semificxible polymers in the melt was determined. The 
stiffness strongly affects the persistence length, end-to-end distance and radius 
of gyration. Stiff chains are more stretched than chains only interacting via 
excluded volume. This also affects the local mutual ordering of the chains. Stiff 
chains pack more parallel on local scales whereas the overall structure remains 
isotropic. The chain length does not influence this strictly local phenomenon. 
Systems consisting of alternating stiff and flexible links behave similar to systems 
with much weaker overall persistence lengths. Their structure is similar to the 
structure of fully flexible chains with larger monomers. Their overall persistence 
length is smaller than expected by a analytical calculation using the persistence 
lengths of the respective potentials. Finally, the overall local structure of the 
melt differs considerably for alternating and homogeneous stiffness. 

The static data presented here shows already ordering effects but its effect 
on dynamical properties relevant to NMR experiments can not be inferred. In 
order to compare directly to NMR experiments on real polymer melts, dynami- 
cal investigations are needed. Such simulations and analyses are presently being 
performed with the mesoscopic model of this article. Moreover, detailed atom- 
istic simulations are underway for melts of specific polymers. They provide 
directly the time evolution of the atom-atom vectors monitored in the experi- 
ments, at least for short times. These simulations will be mapped onto more 
coarse grained simulations like the one presented here. This mapping will al- 
low us to make the connection between the (x, y) parameters of our model and 
real polymers. The presented static properties are a first step on the way of 
understanding ordering phenomena as examined in NMR experiments. 
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